clc
clear all


d = 20; % dimension of permutations
p = 0.5;   % observation probability
q = 0.8; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise
delta = 1;    % noise level in additive noise, delta=0 means no additive noise

Number_of_n = 21;   % number of points
Array_n = 95 + 5 * (1:Number_of_n)';   % number of permutations 

K1 = 5;
K2 = 40;


Number_of_Sim = 30;   % number of repeats

% error 
err_Spectral_K1 = zeros(Number_of_Sim, Number_of_n);
err_QR_K1 = zeros(Number_of_Sim, Number_of_n);
err_Spectral_K2 = zeros(Number_of_Sim, Number_of_n);
err_QR_K2 = zeros(Number_of_Sim, Number_of_n);
err_GPM_K1 = zeros(Number_of_Sim, Number_of_n);
err_GPM_K2 = zeros(Number_of_Sim, Number_of_n);

% running time
time_Spectral_K1 = zeros(Number_of_Sim, Number_of_n);
time_QR_K1 = zeros(Number_of_Sim, Number_of_n);
time_GPM_K1 = zeros(Number_of_Sim, Number_of_n);
time_Spectral_K2 = zeros(Number_of_Sim, Number_of_n);
time_QR_K2 = zeros(Number_of_Sim, Number_of_n);
time_GPM_K2 = zeros(Number_of_Sim, Number_of_n);

for index_n = 1:Number_of_n
    n = Array_n(index_n);    
    
    for index_sim = 1:Number_of_Sim
        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_n, Number_of_n, index_sim, Number_of_Sim);
        iter_time_start = tic;

        [C, row_num, P_true] = Generate_Data_Perm(n, d, delta, p, q);
        
        %% Spectral Estimator
        %[t_Spectral_K1, P_Spectral_K1] = Spectral_Perm_old( C, d);
        
        [t_Spectral_K2, P_Spectral_K2] = Spectral_Perm( C, d, K2);

        %% QR
        %[t_QR_K1, P_QR_K1] = QR_Perm( C, P_Spectral_K1, row_num );

        [t_QR_K2, P_QR_K2] = QR_Perm( C, P_Spectral_K2, row_num );
        
        %% GPM
        %[t_GPM_K1, P_GPM_K1] = GPM_Perm( C, P_Spectral_K1 );

        [t_GPM_K2, P_GPM_K2] = GPM_Perm( C, P_Spectral_K2 );
        
        %%
        fprintf('Iter time: %f\n', toc(iter_time_start));
        %% Update error
        %err_Spectral_K1(index_sim, index_n) = Error_Perm(P_true, P_Spectral_K1);
        %err_QR_K1(index_sim, index_n) = Error_Perm(P_true, P_QR_K1);
        %err_GPM_K1(index_sim, index_n) = Error_Perm(P_true, P_GPM_K1);

        err_Spectral_K2(index_sim, index_n) = Error_Perm(P_true, P_Spectral_K2);
        err_QR_K2(index_sim, index_n) = Error_Perm(P_true, P_QR_K2);
        err_GPM_K2(index_sim, index_n) = Error_Perm(P_true, P_GPM_K2);
        
        %% Update time
        %time_Spectral_K1(index_sim, index_n) = t_Spectral_K1;
        %time_QR_K1(index_sim, index_n) = t_Spectral_K1 + t_QR_K1;
        %time_GPM_K1(index_sim, index_n) = t_GPM_K1 + t_Spectral_K1;

        time_Spectral_K2(index_sim, index_n) = t_Spectral_K2;
        time_QR_K2(index_sim, index_n) = t_Spectral_K2 + t_QR_K2;
        time_GPM_K2(index_sim, index_n) = t_GPM_K2 + t_Spectral_K2;

    end
end

% 
figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_n, 1 - sum(err_Spectral_K2)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_n, 1 - sum(err_GPM_K2)/Number_of_Sim,'-o','LineWidth',2,'MarkerSize',10);
plot(Array_n, 1 - sum(err_QR_K2)/Number_of_Sim,'-*','LineWidth',2,'MarkerSize',10);
%plot(Array_n, 1 - sum(err_QR_K1)/Number_of_Sim,'-*','LineWidth',2,'MarkerSize',10);
%plot(Array_n, 1 - sum(err_Spectral_K1)/Number_of_Sim,'--p','LineWidth',2,'MarkerSize',10);
%plot(Array_n, 1 - sum(err_GPM_K1)/Number_of_Sim,'-^','LineWidth',2,'MarkerSize',10);
hold off
legend('Spectral K=40', 'GPM K=40', 'QR K=40', 'Spectral-K1', 'GPM-K1', 'Location','northwest')
xlabel('\textbf{$n$}','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_n,sum(time_Spectral_K2)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_n,sum(time_GPM_K2)/(Number_of_Sim),'-o','LineWidth',2,'MarkerSize',10);
plot(Array_n,sum(time_QR_K2)/(Number_of_Sim),'-*','LineWidth',2,'MarkerSize',10);
%plot(Array_n,sum(time_QR_K1)/(Number_of_Sim),'-*','LineWidth',2,'MarkerSize',10);
%plot(Array_n,sum(time_Spectral_K1)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
%plot(Array_n,sum(time_GPM_K1)/(Number_of_Sim),'-^','LineWidth',2,'MarkerSize',10);
hold off
legend('Spectral K=40', 'GPM K=40', 'QR K=40', 'Spectral-K1', 'GPM-K1', 'Location','northwest')
xlabel('\textbf{$n$}','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

sgtitle(['$d$ = ', num2str(d), ', $p$ = ', num2str(p), ', $q$ = ', num2str(q), ', $\delta$ = ', num2str(delta)],'Interpreter','latex');

% Save
folder_name = ['k = ', num2str(K2), 'n=', num2str(Array_n(1)),'-',num2str(Array_n(Number_of_n)), '_d=', num2str(d), '_p=', num2str(p), '_q=', num2str(q), '_delta=', num2str(delta)];
mkdir('C:\Users\Administrator\Desktop\Permutation\result',folder_name);
save(['result\',num2str(folder_name),'\err_Spec_K1.mat'], 'err_Spectral_K1');
save(['result\',num2str(folder_name),'\err_LUD_K1.mat'], 'err_QR_K1');
save(['result\',num2str(folder_name),'\err_GPM_K1.mat'], 'err_GPM_K1');
save(['result\',num2str(folder_name),'\err_Spec_K2.mat'], 'err_Spectral_K2');
save(['result\',num2str(folder_name),'\err_LUD_K2.mat'], 'err_QR_K2');
save(['result\',num2str(folder_name),'\err_GPM_K2.mat'], 'err_GPM_K2');
save(['result\',num2str(folder_name),'\time_Spec_K1.mat'], 'time_Spectral_K1');
save(['result\',num2str(folder_name),'\time_LUD_K1.mat'], 'time_QR_K1');
save(['result\',num2str(folder_name),'\time_GPM_K1.mat'], 'time_GPM_K1');
save(['result\',num2str(folder_name),'\time_Spec_K2.mat'], 'time_Spectral_K2');
save(['result\',num2str(folder_name),'\time_LUD_K2.mat'], 'time_QR_K2');
save(['result\',num2str(folder_name),'\time_GPM_K2.mat'], 'time_GPM_K2');
savefig(['result\',num2str(folder_name),'\',num2str(folder_name),'.fig']);